function Fb = bindingModels(xheel,systemData)
% binding model, returns binding force on heel (vector of force in x and y
% dirs) inputs of heel position and the overall 'systemData' structure
% stuff used from systemData
%  systemData.binding.model = string for the binding model ie 'bulldog'
%  systemData.binding.pivotPosition = the current pivot position that the
%    binding is in
%  systemData.binding.preload = length of preload on binding (position
%    difference between natural binding state and when boot is in binding at
%    rest
%  systemData.boot.initialXheel = initial heel position of binding with
%    boot in binding
if strcmpi((systemData.binding.model),'bulldog')
    firstOffset = [64e-3,-23e-3]; %64e-3,-23e-3
    % basic model for now, faking it
    switch systemData.binding.pivotPosition % first is closest to pins and 
        %closest to sole, then moving towards heel, then from pins to heel
        %on the lower row
        case 1
            xp = firstOffset;
        case 2
            xp = firstOffset + 1*[8e-3,0];
        case 3
            xp = firstOffset + 2*[8e-3,0];
        case 4
            xp = firstOffset + [0*8e-3,-8e-3];
        case 5
            xp = firstOffset + [1*8e-3,-8e-3];
        case 6
            xp = firstOffset + [2*8e-3,-8e-3];
        case 7
            xp = firstOffset + [3*8e-3,-8e-3];
        otherwise
            error('Input pivot position not incorporated yet')
    end
    
    % ----------- calc positions --------------------
    
    if isfield(systemData.binding,'preload')
        preload = systemData.binding.preload;
    else
        preload = 0;
    end
    
    xp = xp';
    
    
    initialXheel = systemData.boot.initialXheel;
    
    L0 = norm(initialXheel - xp) - preload; %initial spring assembly length
    % equals the distance between the initial heel position and the pivot,
    % minus whatever preload is being used
    
    Lc = norm(xheel - xp); %current distance between pivot and heel position
    
    hF = [xp-xheel]./norm(xp-xheel); %direction vecotr from heel to pivot
    % gives direction for binding force
    
    dL = Lc - L0; %overall spring length change
    
    
    % ------------- calc spring craziness -------------------------
    
    if length(systemData.binding.springType) ~=length(systemData.binding.springs)
        error('must specify spring lengths and type');
    end
    
    nSprings = length(systemData.binding.springs);
    
    sL = systemData.binding.springs; %lengths of individual spring pieces
    
    for i = 1:nSprings
        % all k values are reference a 4 inch spring
        if strcmpi('red',systemData.binding.springType{i})
            ks(i) = 1.2991e4; %148/2(lb/in)*(kg/2.2lb)(9.81 N/kg)*(in/.0254m)
            bo(i) = .37;
        elseif strcmpi('blue',systemData.binding.springType{i})
            ks(i) = 7.7244e3; %44 lb/in* ....
            bo(i) = 0.5;
        elseif strcmpi('gold',systemData.binding.springType{i})
            ks(i) = 1.8258e4;   %208/2 lb/in *....
            bo(i) = 0.3;
        elseif strcmpi('green',systemData.binding.springType{i})
            ks(i) = 2.5455e4;  %290 lb/in *....
            bo(i) = .25;
        else
            error('Spring type not incorporated yet');
        end
        
        ksPL = ks.*(4*.0254); % spring coefficient times length, sortof nondimensional ish
        
    end
    
    boF = (bo.*sL).*(ksPL./sL); % bottom out force for each spring piece
    
    
    % loop solving spring system and checking for any bottomed out springs
    % repeat until no springs beyond bottom out, then output resulting
    % force
    kBottomOut = 1e6; 
    done = 0;
    cBOS = zeros(1,nSprings); %vector indicating whether that spring piece is bottomed out (1) or not (0)
    
    while ~done
        
        notCBOSInds = find(~cBOS);
        areCBOSInds = find(cBOS);
        
        if length(notCBOSInds) > 1 %if there are more than one spring that isn't bottomed out
            % caculate the effects of each piece, need to solve linear
            % system
            
            % create linear system such that we have something like 
            % the stiffness of each element times its displacement (the
            % unkonwns), are all equal, so do one, minus the next equals
            % zero, then the next minus the one after equals zero and so on
            % for spring pieces that are bottomed out, just take a row, set
            % that spring piece's dx equal to its bottomed out value of dx

            kMat = zeros(nSprings,nSprings);
            fVec = zeros(nSprings,1);

            for i = 1:(length(notCBOSInds)-1)
                kMat(i,notCBOSInds(i)) = ksPL(notCBOSInds(i))./sL(notCBOSInds(i));
                kMat(i,notCBOSInds(i+1)) = -ksPL(notCBOSInds(i+1))./sL(notCBOSInds(i+1));
            end
            
            nInd = length(notCBOSInds)-1;
            for i = 1:(length(areCBOSInds))
                kMat(nInd+i,areCBOSInds(i)) = 1;
                fVec(nInd+i) = bo(areCBOSInds(i))*sL(areCBOSInds(i));
            end
            
            kMat(end,:) = ones(1,nSprings); %last row of system says all dxs must add up to total length change (dL)
            fVec(end) = dL;
            
            dxs = kMat\fVec; %solve linear system to get length changes
            
            newBO = 0;
            for i = 1:length(dxs) %loop over calculated spring length changes
                % check to see if each spring is bottomed out
                
                sF = dxs(i)*ksPL(i)/sL(i); %current spring force
                if sF > boF(i)
                    newBO = 1;
                    cBOS(i) = 1;
                end
            end
            
            if ~newBO % if there were no new bottomed out spring pieces, we are done
                done = 1;
            end
            
            % force magnitude is the same for all the spring pieces, so
            % just take the first not bottomed out spring and get force by
            % multiplying its length change by its stiffness 
            magF = ksPL(notCBOSInds(1))./sL(notCBOSInds(1))*dxs(notCBOSInds(1));
            
        elseif isempty(notCBOSInds) %if there are no springs that are not 
            % bottomed out, find the max bottomed out force from the individual
            % spring pieces, and add that to the bottom out stiffness force
            boDLs = bo(areCBOSInds).*sL(areCBOSInds);
            magF = max(boF) + (dL - sum(boDLs))*kBottomOut;
            done  =1;
        else
            
            done = 1;
            
            if isempty(areCBOSInds) %if there are bottomed out spring pieces
                %add up their bottomed out length changes, if not, just
                %use zero
                boDLs = 0;
            else
                boDLs = bo(areCBOSInds).*sL(areCBOSInds) + 0; 
            end
            
            
            dLi = dL - sum(boDLs); % length change for last spring piece 
            % total length change minus the sum of the length changes from
            % all the bottomed out pieces 
            
%             if dLi < 0
%                 error('Some how bottomed out springs before they should have')
%             end
            Fi = ksPL(notCBOSInds)/sL(notCBOSInds)*dLi; %force on spring from dx
            boL = sL(notCBOSInds)*bo(notCBOSInds); %dx length from bottom out
            
            if dLi > boL %if length change for spring piece is larger than its bottom out length
                % add bottom out stiffness to excess length change
                dxBO = dLi-boL; 
                magF = Fi + dxBO*kBottomOut;
            else
                magF = Fi;
            end
            
        end
                
            
        
        
    end
    
    
   
    Fb = systemData.binding.nTubes*magF.*hF;
    
elseif strcmpi((systemData.binding.model),'hammerhead')
    firstOffset = [-27e-3,-38e-3];
    % basic model for now, faking it
    switch systemData.binding.pivotPosition % first is closest to pins and 
        %closest to sole, then moving towards heel
        case 1
            xp = firstOffset;
        case 2
            xp = firstOffset + 1*[18.8e-3,0];
        case 3
            xp = firstOffset + 2*[18.8e-3,0];
        case 4
            xp = firstOffset + 3*[18.8e-3,0];
        case 5
            xp = firstOffset + 4*[18.8e-3,0];
        otherwise
            error('Input pivot position not incorporated yet')
    end
    
    % ----------- calc positions --------------------
    
    if isfield(systemData.binding,'preload')
        preload = systemData.binding.preload;
    else
        preload = 0;
    end
    
    xp = xp';
    
    
    initialXheel = systemData.boot.initialXheel;
    
    L0 = norm(initialXheel - xp) - preload; %initial spring assembly length
    % equals the distance between the initial heel position and the pivot,
    % minus whatever preload is being used
    
    hF0 = [xp - initialXheel]./norm(xp - initialXheel);
    hF0angle = atan(hF0(2)/hF0(1)); %angle of hF vector off horizontal
    pivotPositionOffset = ...
        (hF0angle < (60*pi/180))*[0.022*sin(hF0angle), 0.022*(1-cos(hF0angle))] + ...
        (hF0angle >= (60*pi/180))*[0.022*sin(60*pi/180), 0.022*(1-cos(60*pi/180))]; 
    offSetCable0 = (hF0angle < (60*pi/180))*.022*hF0angle + (hF0angle >= 60*pi/180)*.022*60*pi/180;
    xp0 = xp + pivotPositionOffset';
    
    L0m = norm(initialXheel - xp0) - preload + offSetCable0;
    
    
    Lc = norm(xheel - xp); %current distance between pivot and heel position
    
    hF = [xp-xheel]./norm(xp-xheel); %direction vector from heel to pivot
    % gives direction for binding force
    
    hFangle = atan(hF(2)/hF(1)); %angle of hF vector off horizontal
    
    if hFangle < (60*pi/180)
        pivotPositionOffset = [0.022*sin(hFangle), 0.022*(1-cos(hFangle))]; %pivot
        %offset due to cable wrap
    elseif hFangle >= (60*pi/180)
        pivotPositionOffset = [0.022*sin(60*pi/180), 0.022*(1-cos(60*pi/180))];
    end
    
    xp = xp + pivotPositionOffset';

    offSetCable = (hFangle < (60*pi/180))*.022*hFangle + (hFangle >= 60*pi/180)*.022*60*pi/180;
    Lcm = norm(xheel - xp) + offSetCable;
    
    hF = [xp-xheel]./norm(xp-xheel); %recalculate direction vector, may want
    % to make this iterative but should be close enough for now.  
   
%     dL = Lc - L0; %overall spring length change
    dL = Lcm - L0m;
    
%     dLt = dL + preload; %overall spring compression 
    
    
    
% ------------- calc spring craziness -------------------------

       if ~isfield(systemData.binding,'springType')
           systemData.binding.springType = 'regular';
       end
    
    
        if strcmpi('regular',systemData.binding.springType)
            ks = 2.7387e4; %156(lb/in)*(kg/2.2lb)(9.81 N/kg)*(in/.0254m)
            bo = 0.048; %bottom out length m
        elseif strcmpi('stiffy',systemData.binding.springType{i})
            ks = 3.423e4; %195 lb/in* .... THIS IS A GUESS NEED TO GET REAL MEASUREMENT
            bo = 0.048; %bottom out length m
        else
            error('Spring type not incorporated yet');
        end
       
    magF = dL*ks + (dL > bo)*(dL - bo)*1e6;

    Fb = magF.*hF;
    
else
    error('Specified binding model not incorporated yet')
end
    